
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%  NOTE THAT THE PRIMARY OUTPUT OF EACH BOOTSTRAP RUN HERE IS JUST STAGE 1 COMMON PARAMETERS 
%%%%%%%  AND A COST FUNCTION.  STUDENT FIXED EFFECT ESTIMATES (PRODUCTIVITY AND MOTIVATION) WILL BE
%%%%%%%  COMPUTED IN PROCESSBOOTSTRAPS.M.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('PointEstimates_MSM7.mat')
UseSeedVectors = 1 ; 
if UseSeedVectors==1 
    load('BootstrapStartValues7.mat') 
end 

BSResults = cell(BootstrapSims,1) ;
%%
diary(strcat('.\BOOTSTRAPS\BSResultsMSM7_',num2str(start),'_',num2str(finish),'_LOG.txt'))




keepnum = 72 ; 
for s=start:finish %BootstrapSims 
    
    timestamp1 = clock ; 
    tic
    disp(strcat('Now computing_',num2str(s-start+1),'_of_',num2str(finish-start+1),'_bootstrap estimates (#',num2str(s),'_of_',num2str(BootstrapSims),')...'))
    datetime(now,'ConvertFrom','datenum') %#ok
    BSResults{s} = struct ;
    sampindx     = [] ;
    for kk=1:length(BSSamp{s,1})
        temp     = find(ProductionDataNonCum.sid==BSSamp{s,1}(kk)) ;  %%%%This is all the rows corresponding to the current re-sampled test subject
        sampindx = [sampindx; [temp, kk*ones(size(temp))]] ; %#ok  %%%%The second column will be a new student id index to replace the old one.  This is necessary since subjects are re-sampled with replacement and therefore can be selected more than once.  Thus, each re-sampled copy of the same student requires a unique new index.
    end
    ProductionSamp     = ProductionDataNonCum(sampindx(:,1),:) ;
    ProductionSamp.sid = sampindx(:,2) ;  %%%%This replaces the old indices with the new indexing for the re-sampled collection of students.
    EPointEsttemp                     = ThetaEEstimator(ProductionSamp,0) ;
    BSResults{s}.EPointEst            = EPointEsttemp ;
    cutoffs                           = EPointEsttemp.cutoffs ; 
        
    sampindx = nan(length(BSSamp{s,1}),2) ;
    for kk=1:length(BSSamp{s,1})
        sampindx(kk,1) = find(AchievementData.id==BSSamp{s,1}(kk)) ;  %%%%This is all the rows corresponding to the current re-sampled test subject
        sampindx(kk,2) = kk ;  %%%%The second column will be a new student id index to replace the old one.  This is necessary since subjects are re-sampled with replacement and therefore can be selected more than once.  Thus, each re-sampled copy of the same student requires a unique new index.
    end
    AchievementSamp           = AchievementData(sampindx(:,1),:) ;
    AchievementSamp.id        = sampindx(:,2) ;  %%%%This replaces the old indices with the new indexing for the re-sampled collection of students.
    AchievementSamp.ThetaE(:) = NaN ;
    for kk=1:length(BSResults{s}.EPointEst.studlistFULL)
        AchievementSamp.ThetaE(AchievementSamp.id==BSResults{s}.EPointEst.studlistFULL(kk)) = BSResults{s}.EPointEst.ThetaE(kk) ;
    end

    knotvec_c = LPointEst.knotvec_c;
    P = ThetaLEstimatorMSM(EPointEsttemp,AchievementSamp,knotvec_c,[],LPointEst.params,0,...
        [1;max(LPointEst.Q80freq(:,1,1));max(LPointEst.Q80freq(:,1,2));max(LPointEst.Q80freq(:,1,3))],1) ;
    
    BSstate      = [12,max(LPointEst.QuantileFunq1),max(LPointEst.QuantileFunq2),max(LPointEst.QuantileFunq3)] ;
    
    pi_c = LPointEst.params ; 
    if UseSeedVectors==1 
        pigrid = nan(length(Solutions{1,1}),length(Solutions)) ; 
        for ii=1:length(Solutions) 
            pigrid(:,ii) = Solutions{ii,1} ; 
        end 

        pigrid = [pi_c pigrid] ; %#ok 
    else
        pigrid = pi_c ;
    end
    
    [pigrid,ssrlog,survivorfrac] = ASPS(pigrid,knotvec_c,P,BSstate,keepnum,1) ;
    survivorfrac  %#ok
    [~,indx] = mink(ssrlog,1) ;
    pi_c0 = pigrid(:,indx) ;


    disp('%%%%PSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPS%%%%%')
    disp(strcat('%%%% PRELIMINARY ROUND OF SIMULATED GMM ESTIMATION VIA PATTERN SEARCH...     %%%%%'))
    disp('%%%%PSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPS%%%%%')
    BSstate(1) = 2 ;
    LPointEsttemp = ThetaLEstimatorMSM(EPointEsttemp,AchievementSamp,knotvec_c,[],pi_c0,0,BSstate,0,pigrid') ;
    try  %%%%Once in a while, the first round of pattern search returns something infeasible, in which case ASPS will
        %%%% throw an error.  So, for that we have this try-catch statement to keep the process moving.
        [pigrid,ssrlog,survivorfrac] = ASPS([LPointEsttemp.params pigrid],knotvec_c,P,BSstate,keepnum,0) ;
        survivorfrac  %#ok
    catch
        disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
        disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PATTERN SEARCH RETURNED SOMETHING INFEASIBLE, SO WE RUN A PRELIMINARY ROUND OF GA !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
        disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
        BSstate(1) = 12 ;
        disp('%%%%GAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGA%%%%%')
        disp(strcat('%%%% SIMULATED GMM ESTIMATION, PARTIAL ROUND VIA GENETIC ALGORITHM...'))
        disp('%%%%GAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGA%%%%%')
        LPointEsttemp = ThetaLEstimatorMSM(EPointEsttemp,AchievementSamp,knotvec_c,[],[],0,BSstate,0,pigrid') ;

        [pigrid,ssrlog,survivorfrac] = ASPS([LPointEsttemp.params; LPointEsttemp.population]',knotvec_c,P,BSstate,keepnum,0) ;
        survivorfrac  %#ok
    end


    BSstate(1) = 13 ;
    disp('%%%%GAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGA%%%%%')
    disp('%%%%  TERMINAL ROUND OF SIMULATED GMM ESTIMATION VIA GENETIC ALGORITHM...  %%%%%')
    disp('%%%%GAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGA%%%%%')
    LPointEsttemp = ThetaLEstimatorMSM(EPointEsttemp,AchievementSamp,knotvec_c,[],[],0,BSstate,0,pigrid') ;
    ssrGAterminal = LPointEsttemp.ssrvalc ;
    pi_c0GA = LPointEsttemp.params' ;

    
    
    disp('%%%%PSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPS%%%%%')
    disp(strcat('%%%% TERMINAL ROUND OF SIMULATED GMM ESTIMATION VIA PATTERN SEARCH...     %%%%%'))
    disp('%%%%PSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPS%%%%%')
    BSstate(1) = 3 ;
    LPointEsttemp = ThetaLEstimatorMSM(EPointEsttemp,AchievementSamp,knotvec_c,[],LPointEsttemp.params',0,BSstate,0,pigrid') ;
    ssrPSterminal = LPointEsttemp.ssrvalc ;
    format long
    LPointEsttemp.params
    BSResults{s}.LPointEst = LPointEsttemp ;

    BSResults{s}.LPointEst.ssrPSpenultimate = ssrPSterminal ;
    timestamp2             = clock ;
    ElapsedMinutes         = etime(timestamp2,timestamp1)/60 %#ok

    if mod(s,10)==0  %%%% After every 10 units, we pause to save progress...
        disp('***')
        disp('***')
        disp('%%%%%%%%%%%%%%%%%%%%SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE%%%%%%%%%%%%%%%%%%%%')
        disp('%%%%%%%%%%%%%%%%%%%%SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE%%%%%%%%%%%%%%%%%%%%')
        disp('%%%%%%%%%%%%%%%%%%%%SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE%%%%%%%%%%%%%%%%%%%%')
        disp('%%%%%%%%%%%%%%%%%%%%                                                                                                                       %%%%%%%%%%%%%%%%%%%%')
        disp('%%%%%%%%%%%%%%%%%%%%      PAUSING WORK BRIEFLY TO SAVE PROGRESS, SO THE IT JERKS DON''T KEEP SCREWING US WITH RANDOM RE-BOOTS...           %%%%%%%%%%%%%%%%%%%%')
        disp('%%%%%%%%%%%%%%%%%%%%                                                                                                                       %%%%%%%%%%%%%%%%%%%%')
        disp('%%%%%%%%%%%%%%%%%%%%SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE%%%%%%%%%%%%%%%%%%%%')
        disp('%%%%%%%%%%%%%%%%%%%%SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE%%%%%%%%%%%%%%%%%%%%')
        disp('%%%%%%%%%%%%%%%%%%%%SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE_SAVE%%%%%%%%%%%%%%%%%%%%')
        disp('***')
        disp('***')
        command = strcat('BSResultsMSM_',num2str(start),'_',num2str(finish),'=BSResults ;') ;
        eval(command) ;
        save(strcat('.\BOOTSTRAPS\BSResultsMSM7_',num2str(start),'_',num2str(finish),'.mat'),strcat('BSResultsMSM_',num2str(start),'_',num2str(finish)),'-v7.3')
    end
end






toc
diary off



